//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: alexander.howard@cern.ch
// --------------------------------------------------------------
// Comments
//
//                  Underground Advanced
//               by A. Howard and H. Araujo 
//                    (27th November 2001)
//
//
// PrimaryGeneratorAction program
// --------------------------------------------------------------

#include "DMXPrimaryGeneratorAction.hh"
#include "G4ParticleGun.hh"
#include "G4Event.hh"
#include "G4Run.hh"
#include "G4UnitsTable.hh"
#include "G4ParticleTable.hh"
#include "Randomize.hh"
#include "G4RandomDirection.hh"

#include "globals.hh"


DMXPrimaryGeneratorAction::DMXPrimaryGeneratorAction(G4int RunID) 
  : G4VUserPrimaryGeneratorAction(),particleGun(0)
{
  runid = RunID;
  G4int Nparticle = 1;  
  particleGun = new G4ParticleGun(Nparticle);
}
DMXPrimaryGeneratorAction::DMXPrimaryGeneratorAction() 
  : G4VUserPrimaryGeneratorAction(),particleGun(0)
{
  runid = 0 ;
  G4int Nparticle = 1;  
  particleGun = new G4ParticleGun(Nparticle);
}


DMXPrimaryGeneratorAction::~DMXPrimaryGeneratorAction() {

  delete particleGun;
  delete particle;
  delete particleTable;
}


void DMXPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {

  G4int eventid = anEvent->GetEventID();
  G4int iz   = runid / 305;

  G4int nn;
  //  if( iz<24 ) nn=iz*37500;
  //  else        nn=900000;

  //nn=900000;
  //G4int NOP = (100000+nn)/10;  //Set Number Of Particle per 1 of 10Event
  G4int NOP = 1;

  for(int i=0;i<NOP;i++){

  //--- Ion
  //SetParticleIon( eventid );
  //--- Photon
  //  SetParticlePhoton(runid,i);
  //--- Neutron
  //SetParticleNeutron(eventid);
  //--- Gamma
  SetParticleGamma(eventid);
  //--- Cosmic
  //  SetParticleCosmic();
    //    if(i%10000==0)    G4cout<<"ParticleGun:"<<i<<G4endl;
    particleGun->GeneratePrimaryVertex( anEvent );
  }

  energy_pri = 0.;
  energy_pri = particleGun->GetParticleEnergy();

}

void DMXPrimaryGeneratorAction::SetParticleIon(G4int ievt){

  //---Ion
  //  G4String particleName;
  //  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  //alpha---
  //particle = particleTable->FindParticle(particleName="alpha");
  //particleGun->SetParticleDefinition(particle);

  //Ar ion---
  G4int iZ = 18;
  G4int iA = 40;
  G4float iC = 1.*eplus;
  G4float excitE = 0.*keV;
  particle = particleTable->GetIon(iZ,iA,excitE);
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleCharge(iC);
  
  //----- Set Energy
  particleGun->SetParticleEnergy(300.*keV);

  //----- Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  G4double phi  = 2*3.141592*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4float xx=sint*cosp;
  G4float yy=sint*sinp;
  G4float zz=cost;
  G4ThreeVector vecgun(xx,yy,zz);
  //G4ThreeVector vecgun(0.,0.,-1.);
  particleGun->SetParticleMomentumDirection( vecgun );

  ievt = ievt / 100;
  G4int iz = ievt / 21 % 18;
  G4int ir = ievt % 21;

  G4ThreeVector vecpos(-120.+ir*10.*mm, 0., -545.+iz*10.*mm);
  //  G4ThreeVector vecpos(0., 0., -505.*mm);
  particleGun->SetParticlePosition( vecpos );

  //G4cout<<"Particleposition("<<-120.+ir*10.*mm<<","<<-545.+iz*10.*mm<<")"<<G4endl;

}

void DMXPrimaryGeneratorAction::SetParticlePhoton(G4int ievt,G4int iphoton){

  //---photon
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  
  //---Set Photon Energy
  G4float PhotonEnergy = 9.69*eV;
  PhotonEnergy = GetPhotonEnergy();

  particle = particleTable->FindParticle(particleName="opticalphoton");
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleEnergy(PhotonEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  G4double phi  = 2*3.141592*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  particleGun->SetParticleMomentumDirection( vecgun );

  //---Set Position
  /*
  G4int ix = ievt % 30;
  G4int iy = ievt / 30 % 30;
  G4int iz = ievt / 30 / 30 % 34;
  G4int xx = (-145.+ix*10.)*mm;
  G4int yy = (-145.+iy*10.)*mm;
  G4int zz;
  if( iz<24 )  zz = (-595.+iz*10.)*mm;
  else         zz = (-360.+(iz-24)*1.)*mm;
  */
  /*
  G4int ir   = ievt % 15;
  G4int iphi = ievt / 15 % 13;
  G4int iz = ievt / 15 / 13 % 34;
  G4float rr    = (5.+ir*10.)*mm;
  G4float pphi  = iphi*15.*deg;
  G4float zz;
  if( iz<24 )  zz = (-595.+iz*10.)*mm;
  else         zz = (-360.+(iz-24)*1.)*mm;
  */


  int nphi[15] = {3,5,8,11,14,17,20,24,27,30,33,36,39,19,19}; //305 
  int id=-1;

  int runID=ievt % 305;
  if( runID>=0  && runID<nphi[0] )  id=0;
  else if( runID<(nphi[0]+nphi[1]) )     id=1;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]) )  id=2;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]) )  id=3;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]) )  id=4;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]) )  id=5;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]) )  id=6;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]) )  id=7;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]) )  id=8;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]) )  id=9;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]) )  id=10;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]) )  id=11;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]) )  id=12;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]+nphi[13]) )  id=13;
  else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]+nphi[13]+nphi[14]) )  id=14;
  else if( id==-1 ){ G4cout<<runid<<G4endl;return;}

  int sum=0;
  if(id!=0) for(int ii=0;ii<id;ii++) sum+=nphi[ii];
  int iphi = (runID-sum) % nphi[id];
  int ir   = id;
  int iz   = ievt / 305 ;

  G4float rr    = (5.+ir*10.)*mm;
  G4float pphi  = iphi*(180/(G4float)(nphi[ir]-1))*deg;
  G4float zz;
  //  if( iz<24 )  zz = (-595.+iz*10.)*mm;
  //  else         zz = (-360.+(iz-24)*1.)*mm;
  zz = (-360.+iz*1.)*mm;


  G4ThreeVector vecpos( rr*cos(pphi), rr*sin(pphi), zz);
  //  G4ThreeVector vecpos( xx, yy, zz);

  particleGun->SetParticlePosition( vecpos );
  if(iphoton==0){
    //G4cout<<"Photon Energy = "<<G4BestUnit(PhotonEnergy,"Energy")<<G4endl;
    //    G4cout<<"Position("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
    //    G4cout<<ir<<" "<<iphi<<" "<<rr<<" "<<pphi*180/3.14<<G4endl;
    G4cout<<"Position("<<rr*cos(pphi)<<","<<rr*sin(pphi)<<","<<zz<<")"<<G4endl;
  }
}

void DMXPrimaryGeneratorAction::SetParticleNeutron(G4int ievt){

  //---Neutron
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="neutron");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float NeutronEnergy;
  //NeutronEnergy = Get241AmBeEnergy()*MeV; // 241Am-Be(Bare)
  //NeutronEnergy = Get252CfEnergy()*MeV;   // 252Cf

  NeutronEnergy = 1.*MeV;
  particleGun->SetParticleEnergy(NeutronEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  G4double phi  = 2*3.141592*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  //  G4ThreeVector vecgun(1.,0.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  //---Set Position
  G4int iz = ievt / 1000;
  //  G4float xx = -172.*mm;
  G4float xx = -180.*mm;
  G4float yy = 0.*mm;
  G4float zz = (-545.+iz*20.)*mm;

  //G4ThreeVector vecpos(xx, yy, zz);
  //  G4cout<<"Particleposition("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  G4cout<<"NeutronEnergy = "<<G4BestUnit(NeutronEnergy,"Energy")<<G4endl;
  G4ThreeVector vecpos(xx, yy, -455.*mm);

  particleGun->SetParticlePosition( vecpos );
}

void DMXPrimaryGeneratorAction::SetParticleGamma(G4int ievt){
  //---Gamma
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="gamma");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float GammaEnergy;
  GammaEnergy = 1.33*MeV;
  particleGun->SetParticleEnergy(GammaEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  G4double phi  = 3.14159265358979*(G4UniformRand()+1/4.);
  //  G4double phi  = 0.5*3.141592*(1.-2.*G4UniformRand());
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  // G4ThreeVector vecgun(-1.,1.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  //---Set Position
  G4int iz = ievt / 10000;
  G4float PHI = -67.5*deg;
  G4float xx = 180.*cos(PHI)*mm;
  G4float yy = 180.*sin(PHI)*mm;
  G4float zz = -465.*mm;
  //G4float zz = (-545.+iz*20.)*mm;

  //G4ThreeVector vecpos(xx, yy, zz);
  //  G4cout<<"Particle Position("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  G4cout<<"      Beam Vector("<<sint*cosp<<","<<sint*sinp<<","<<cost<<")"<<G4endl;
  //  G4cout<<"Gamma Energy = "<<G4BestUnit(GammaEnergy,"Energy")<<G4endl;
  G4ThreeVector vecpos(xx, yy, zz);

  particleGun->SetParticlePosition( vecpos );

}

void DMXPrimaryGeneratorAction::SetParticleCosmic(){
  //---Cosmic
  gRandom->SetSeed(time(NULL));
  TFile *fread=new TFile("./inputFiles/MuonFlux.root");
  TH1F *hSpectrumMu_p=(TH1F*)fread->Get("hSpectrumMu_p");
  TH1F *hSpectrumMu_n=(TH1F*)fread->Get("hSpectrumMu_n");
  TH1F *hSpectrumMu  =(TH1F*)fread->Get("hSpectrumMu");

  RatioPtoALL = hSpectrumMu_p->Integral()/hSpectrumMu->Integral();

  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  float vRandom = gRandom->Rndm();
  float MuonEnergy = 0;

  if( RatioPtoALL>vRandom ){
    particle = particleTable->FindParticle(particleName="mu+");
    MuonEnergy = hSpectrumMu_p->GetRandom();
  }
  else{
    particle = particleTable->FindParticle(particleName="mu-");
    MuonEnergy = hSpectrumMu_n->GetRandom();
  }

  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleEnergy(MuonEnergy*MeV);

  G4cout<<"Ratio = "<<RatioPtoALL<<G4endl;
  G4cout<<"MuonType   = "<<particleName<<G4endl;
  G4cout<<"MuonEnergy = "<<G4BestUnit(MuonEnergy,"Energy")<<G4endl;

  theta = G4UniformRand()*3.14;
  phi   = (G4UniformRand()-0.5)*3.14;

  G4ThreeVector vecgun(-1.,0.,0.);
  G4ThreeVector vecpos(250.*mm,0.*mm,-500.*mm);
  //  G4ThreeVector vecgun(-1.,0.,-0.5);
  //G4ThreeVector vecpos(250.*mm,0.*mm,200.*mm);
  particleGun->SetParticleMomentumDirection( vecgun );
  particleGun->SetParticlePosition( vecpos );

  fread->Delete();
  hSpectrumMu_p->Delete();
  hSpectrumMu_n->Delete();
  hSpectrumMu  ->Delete();
}

G4double DMXPrimaryGeneratorAction::GetPhotonEnergy(){
  TF1 *fAr=new TF1("fAr","exp(-0.5*((x-128.0)/2.929)*((x-128.0)/2.929))",0,400.);
  G4double PlankH = 4.135667516e-15*eV*s;
  G4double LightC = 2.99792458e8*m/s; 
  G4double PhotonEnergy = PlankH*LightC/(fAr->GetRandom()*nm);
  fAr->Delete();

  return PhotonEnergy;
}
G4double DMXPrimaryGeneratorAction::Get252CfEnergy(){
  TF1 *fNspectrum=new TF1("fNspectrum","sqrt(x)*exp(-x/1.42)",0.01,100);
  G4double NeutronEnergy = fNspectrum->GetRandom();
  fNspectrum->Delete();

  return NeutronEnergy;
}
G4double DMXPrimaryGeneratorAction::Get241AmBeEnergy(){
  TH1F *h1=new TH1F("h1","h1",10000,0,100);

  ifstream fin("./inputFiles/AmBeTable.dat");
  if(!fin){G4cout<<"No Such File"<<G4endl;return 0;}

  while( !fin.eof() ){
    float Energy,Cf,AmBe,tmp;
    fin >> Energy >> Cf >> tmp >> AmBe >> tmp;
    h1->Fill((Energy/1.0e6),AmBe);
  }
  int nbin = h1->GetNbinsX();
  for(int ibin=0;ibin<nbin;ibin++){
    float vbin = h1->GetBinContent(ibin+1);
    if( ibin<15.8/(100/10000.)&&!vbin&&ibin>0 ){
      h1->SetBinContent(ibin+1,h1->GetBinContent(ibin));
    }
  }
  G4double NeutronEnergy = h1->GetRandom();
  h1->Delete();

  return NeutronEnergy;
}
